cd "/Users/delgerjargaluvsh/Dropbox/Synthesis project docs/russia_forest/Code and Data/Uploaded to Dataverse" 
clear	

version 13.1
***********************************************************************************
***********************************************************************************
************Merging satellite data into one master dataset 
*Load dynamic forest types data
use FCC_satellitedata_dynamic_types_ha.dta 
*Merge with forest gain data 
merge 1:1 newid using FCC_satellitedata_forest_gain_ha.dta
drop _merge 
*Merge with forest loss data 
merge 1:1 newid using FCC_satellitedata_loss_type_alloc_ha.dta
drop _merge 
*Drop empty columns
drop Q
drop R 

*Merge with forest cover data 
merge 1:1 newid using stat_with_forest_cover_1.dta
drop _merge

***********Dropping observations that are not in Russia from the master dataset 
keep if newid < 5 | newid > 75 
drop if newid > 132

*Merge the Soil index data with the main dataset 
merge 1:1 newid using Soil_index_data.dta
drop _merge 

*Merge the Evergreen forest data with the main dataset 
merge 1:1 newid using final_forest_types_Rus_regions_data.dta
drop _merge 

*Saving the data as mastersheet 
save "mastersheet"

************Adding regional income data 
********Loading income data for 2001 (income_1991_2004.dta) 
clear
use income_1991_2004.dta 
*Keep data for year 2001 (89 regions)
keep if year == 11

*Merge regional nominal income data with CPI adjusted income data (used later for robustness check)
merge 1:1 id using income_nominal_adjustment_2001.dta
drop _merge 

*Drop observations sharing ter_ao identifier with others 
*Evenkiisky; Krasnoyarsky; Taimirskiy
drop if ter_ao == 1104 
*Ust'-Ordinskiy Buryatskiy; Irkutskaya oblast 
drop if ter_ao == 1125
*Kamchatskaya oblast; Koryaksky a.o.
drop if ter_ao == 1130
*Aginskiy Buryatskiy a.o.; Chitiskaya obl.
drop if ter_ao == 1176
*Drop Komi-Permyak, because regional income value for Perm' includes that of Komi-Permyak 
sort ter_ao
drop if id==81 
*Drop Nenets, because regional income value for Arkhangel'sk includes that of Nenets 
drop if id==83
*Drop Tyumen, which is outside the study area 
drop if id==72 
*Khanty-Mansk and Yamalo-Nenets have the same ter_ao. So, replace ter_ao for Khanty-Mansk with 9999. 
replace ter_ao=9999 in 51
sort ter_ao
save "income2001" 
*Regional income data for 2001 with unique ter_ao is saved as income2001

*********Loading income data for 1990
clear
use income_1985_1998.dta
*Keep data for year 1990 
keep if year== 90
*Rename ter (territory id) variable as ter_ao 
rename ter ter_ao
*Merge the regional income data for 1990 with regional income data for 2001 
merge 1:1 ter_ao using income2001.dta
rename _merge merge_income90and01

**********Merging income data with sheet with newid variable and region identifiers 
merge 1:1 ter_ao using FCCproject_identifiers.dta
rename _merge merge_income90and01andidentifier
save "income1990_2001_identifiers" 

**********Loading income data for 1991 for Adygey and Karachay Cherkess 
clear
use income_1985_1998.dta
*Keep data for year 1991 for Adygey and Karachay Cherkess
keep if year== 91
keep if ter==1179|ter==1191
*Rename ter (territory id variable) as ter_ao 
rename ter ter_ao
merge 1:1 ter_ao using income1990_2001_identifiers.dta

*Rename the regional income pc variables 
rename mincomepop_percap living_st2001
rename nincpc living_st1990
rename adjusted_income adjusted_living_st2001

*Drop regions not included in the study area 
drop if newid >= . 
drop _merge

*Merge the regional income pc datasheet with the master datasheet 
merge 1:1 newid using mastersheet.dta
drop _merge 

*Merge the master datasheet with dataset of additional covariates (distance nearest city, governance scores, autonomous status, % urban population, access, bureaucrats)
merge 1:1 newid using FCCproject_IV.dta
drop _merge 
drop id id_ao year 

************Combining the data for Komi-Permyak and Perm' as well as Arkhangelsk and Nenets************
****Modify the data on Perm' to include Komi-Permyak 
**Collapse Komi-Permyak into Perm' 
*Change newid of Komi-Permyak with the same newid as Perm' (Both Perm' (line/obs 43) and Komi-Permyak (line/obs 25) have newid = 114)
replace newid=newid[43] in 25
*Command to sum/collapse all variables by newid (Instead of 61, now we have 60 obs)
*Bring variables in string format to the beginning of the dataset 
order newid id_code nameeng merge_income90and01 Region id_code merge_income90and01andidentifier cntry S ID_CODE CNTRY approved_names Country Province officials_1995 officials_1996 officials_1997 officials_1998 officials_1999 workers_1995 workers_1996 workers_1997 workers_1998 workers_1999 , first
collapse (sum) ter_ao-bur_2011, by(newid) 
****Variables from the satellite data are all summed on the level of Perm'. 
*The variables that are summed and are subsequently not altered include forest dynamic types, forest gain and forest loss due to other and fire. 
*Other summed variables include per_data and forest_ha, but Perm' and Komi_Permyak are both 100% covered and have more than 10k ha of forest. 
*Some covariates are also summed over the two regions and are not altered, including area under soil indices (water, all_classes), forest bureaucrats and total forest area. 
*Regional income per capita are values for Perm', because for 1990/1991 regional income pc is reported only for Perm'. The 2001 figure for Perm' includes Komi Permyak. 
*Autonomous status variable from Perm' is used.  
*However, some covariates need to be adjusted as below. 
*Accessibility variable is recalculated in ArcGIS for larger Perm' and entered into the datasheet with covariates for both regions. 
*Since values for both regions are collapsed, the value for Komi Permyak is deducted. 
replace avv50k_mn = avv50k_mn[42]-340 in 42
*Recalculate the percentage evergreen forest variable 
replace pr_evergreen = EvergreenNeedleleafForest[42]/totalforestedarea[42] in 42
*In order to use distance nearest city value for Perm' for this observation, we deduct the value of Komi-Permyak
replace dist_nearest_city = dist_nearest_city[42]-1394 in 42
*In order to calculate pr_upop value for Perm for this observation, we deduct the value of Komi-Permyak, because values from the Stat yearbooks include Komi-Permyak. 
replace pr_upop1991 = pr_upop1991[42]-30 in 42
replace pr_upop2001 = pr_upop2001[42]-26.1 in 42
*Average governance scores (Perm' has score of 41 for both periods; Komi Permyak has scores of 26 and 29 for the 1990s and 2000s respectively)
replace governance_91_01 = governance_91_01[42]/2 in 42
replace governance = governance[42]/2 in 42
*Changing autonomous status variable to Perm's value. 
replace autonom = 0 in 42

****Modifying the data on Arkhangel'sk to include Nenets************
**Collapsing Nenets into Arkhangel'sk 
*Change newid of Nenets with the same newid as Arkhangel'sk (Both Arkhangel'sk (line/obs 8) and Nenets (line/obs 35) have newid = 79)
replace newid=newid[8] in 35
*Command to sum/collapse all variables by newid (Instead of 60, now we have 59 obs)
order newid, first
collapse (sum) ter_ao-bur_2011, by(newid) 
****Variables from the satellite data are all summed on the level of Arkhangelsk. 
*The variables that are summed and are subsequently not altered include forest dynamic types, forest gain and forest loss due to other (logging) and fire. 
*Other summed variables include per_data and forest_ha, but Arkhangel'sk and Nenets are both 100% covered and have more than 10k ha of forest. 
*Some covariates are also summed over the two regions and are not altered, including area under soil indices (water, all_classes), forest bureaucrats and total forest area. 
*Regional income per capita are values for Arkhangel'sk, because for 1990/1991 regional income pc is reported only for Arkhangel'sk. The 2001 figure for Arkhangel'sk includes Nenets. 
*Autonomous status variable from Arkhangel'sk is used.  
*However, some covariates need to be averaged as done below.
*Accessibility variable is recalculated in ArcGIS for the larger region Arkhangel'sk and entered into the datasheet with covariates for both regions. Since values for both regions are collapsed, the value for Nenets is deducted. 
replace avv50k_mn = avv50k_mn[8]-1073 in 8
*Recalculate the percentage evergreen forest variable 
replace pr_evergreen = EvergreenNeedleleafForest[8]/totalforestedarea[8] in 8
*In order to use distance nearest city value for Arkhangel'sk for this observation, we deduct the value of Nenets. 
replace dist_nearest_city = dist_nearest_city[8]-2230 in 8
*In order to include pr_upop value for Arkhangel'sk for this observation, we deduct the value of Nenets. Values from the Stat yearbooks include Nenets. 
replace pr_upop1991 = pr_upop1991[8]-63.1 in 8
replace pr_upop2001 = pr_upop2001[8]-60.4 in 8
*Average governance scores (Arkhangel'sk has score of 37 for both periods; Nenets has scores of 23 and 27 for 1990s and 2000s respectively)
replace governance_91_01 = governance_91_01[8]/2 in 8
replace governance = governance[8]/2 in 8
*Changing autonomous status variable to Perm's values 
replace autonom = 0 in 8

*****************************************************************************
************Generating the dependent variables************
****Figures for normalization 
*1985 total forest (stable forest + forest loss + forest loss followed by forest gain + repeated forest loss separated by forest gain)
gen totalforest85 = dyn_typ_b + dyn_typ_d + dyn_typ_e + dyn_typ_f
*2000 total forest (stable forest + forest gain over non-forest in 1985 + repeated forest loss separated by forest gain + forest loss on areas which gain forest cover after non-forest state in 1985)
gen totalforest00 = dyn_typ_b + dyn_typ_c + dyn_typ_f + dyn_typ_g
*Alternative measure of 2000 total forest (totalforest85 + 2000 gains - 2000 loss) is forest_cover_2000

*Logging loss
gen loggingloss_1990s = otherloss_1989 + otherloss_1990 + otherloss_1991 + otherloss_1992	+ otherloss_1993 + otherloss_1994 + otherloss_1995 + otherloss_1996	+ otherloss_1997 + otherloss_1998 + otherloss_1999 + otherloss_2000
gen loggingloss_2000s = otherloss_2001 + otherloss_2002 + otherloss_2003 + otherloss_2004	+ otherloss_2005 + otherloss_2006 + otherloss_2007 + otherloss_2008	+ otherloss_2009 + otherloss_2010 + otherloss_2011 + otherloss_2012

gen pr_logging_1990s = (loggingloss_1990s/totalforest85)*100
gen pr_logging_2000s = (loggingloss_2000s/totalforest00)*100
gen alt_pr_logging_2000s = (loggingloss_2000s/forest_cover_2000)*100

* Fire loss 
gen fireloss_1990s = fireloss_1989 + fireloss_1990 + fireloss_1991 + fireloss_1992 + fireloss_1993 +	fireloss_1994 + fireloss_1995 + fireloss_1996 + fireloss_1997 + fireloss_1998 + fireloss_1999 + fireloss_2000
gen fireloss_2000s = fireloss_2001 + fireloss_2002 + fireloss_2003	+ fireloss_2004	+ fireloss_2005 +	fireloss_2006	+ fireloss_2007	+ fireloss_2008	+ fireloss_2009	+ fireloss_2010	+ fireloss_2011	+ fireloss_2012

gen pr_fireloss_1990s = (fireloss_1990s/totalforest85)*100
gen pr_fireloss_2000s = (fireloss_2000s/totalforest00)*100
gen alt_pr_fireloss_2000s = (fireloss_2000s/forest_cover_2000)*100

*Gain over cropland or pasture
gen totalarea_km = (all_classes - water)
*Check it against the GIS area_km2 data 
corr totalarea_km area_km2
*Convert to ha 
gen totalarea = totalarea_km*100 

gen total_nonforest85 = totalarea - totalforest85
gen total_nonforest00 = totalarea - totalforest00
gen alt_total_nonforest00 = totalarea - forest_cover_2000

gen pr_croppastgain1990s =(((aff_2000_crop_pas/16)*12)/total_nonforest85)*100
gen pr_croppastgain2000s = (aff_2012_crop_past/total_nonforest00)*100
gen alt_pr_croppastgain2000s = (aff_2012_crop_past/alt_total_nonforest00)*100

************Generating covariates************
*Scale distance to nearest major city per 1000 km 
gen rescaled_dist = dist_nearest_city/1000
*Calculate fertile area (crop suitability) as percentage of all land, except lakes/rivers
gen fertile_area = ((si41_55 + si56_70 + si71_85 + si86_100)/(all_classes - water))*100
*Take log of regional income per capita variables 
gen ln_lst1990 = log(living_st1990)
gen ln_lst2001 = log(living_st2001)
gen ln_adj_lst2001 = log(adjusted_living_st2001) 
*Take log of total forest variables 
gen ln_totalfor85 = log(totalforest85)
gen ln_totalfor00 = log(totalforest00)
gen alt_ln_totalfor00 = log(forest_cover_2000)

*Convert remoteness to nearest city with at least 50k people from minutes to hours 
gen access = (avv50k_mn/60)

*Take squares of governance variables 
gen gov_sq_90 = governance_91_01*governance_91_01
gen gov_sq_00 = governance*governance

*Scale percent evergreen forest variable 
gen pr_evergreen_for = pr_evergreen*100

****Normalize and take log of forest bureaucrats variables 
*1988 all bureaucrats (For Adegeye, we use 1992 bureaucrats data)
gen norm_off1988 = allworkers_1988/totalforest85
gen ln_off1988 = log(norm_off1988)

*1992 all bureaucrats 
gen norm_off1992 = allworkers_1992/totalforest85
gen ln_off1992 = log(norm_off1992)

*2001 all bureaucrats
gen norm_off2001 = allworkers_2001/totalforest00
gen ln_off2001 = log(norm_off2001)

*Normalize with alternative total forest figures 
gen alt_norm_off2001 = allworkers_2001/forest_cover_2000
gen alt_ln_off2001 = log(alt_norm_off2001)


************Checking for missingness in the data************
*Regions with missing/unreported 1990 regional income pc values include Chechnya, Komi-Permyak, Khanty-Mansk, Nenets, Yamalo-Nenets (For Adegey and Karachay Cherkess, we use 1991 values)
*Regions with missing/unreported 2001 regional income pc values include Chechnya. 
**Regions with missing/unreported 1988 forest bureaucrats values include Chechnya, Komi-Permyak, Nenets, Khanty-Mansk, Yamalo-Nenets, and Moscow (For Adegey, we use 1992 bureaucrats data)
**Regions with missing/unreported 2001 forest bureaucrats values include Chechnya, Ingush, Yamalo-Nenets, Nenets, Moscow, and Komi-Permyak. 
*** Chechnya is missing governance score variables. 

mdesc pr_logging_1990s pr_logging_2000s pr_fireloss_1990s pr_fireloss_2000s pr_croppastgain1990s pr_croppastgain2000s fertile_area rescaled_dist pr_evergreen_for access pr_upop2001 pr_upop1991 ln_lst1990 ln_lst2001 autonom ln_off1988 ln_off2001 governance_91_01 governance

************Creating the sample and summary statistics table************ 
*Drop the regions that do not have enough data coverage (Khanty-Mansk, Yamalo-Nenets, Sverdlovsk, and Chelyabinsk)
drop if per_data < 100

*Drop the regions with less than 10k hectares of forest (Astrakhan, Kalmykia)
drop if forest_ha < 10000

* Drop St. Petersburg
drop if newid == 86 

*Create the sample (49 obs)
regress pr_logging_1990s pr_logging_2000s pr_fireloss_1990s pr_fireloss_2000s pr_croppastgain1990s pr_croppastgain2000s fertile_area rescaled_dist pr_evergreen_for access pr_upop1991 pr_upop2001 ln_lst1990 ln_lst2001 autonom ln_off1988 ln_off2001 governance_91_01 governance
estimates store m2
generate sample = e(sample)
drop if sample == 0

merge 1:1 newid using FCCproject_identifiers.dta
drop _merge 
save "key_variables" 


*Summary statistics (Table 1)
*Table 1 
sutex fertile_area access pr_evergreen_for rescaled_dist autonom pr_logging_1990s pr_fireloss_1990s pr_croppastgain1990s pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 pr_logging_2000s  pr_fireloss_2000s pr_croppastgain2000s pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance, minmax 

*************************************************************
*************************************************************
***********Connecting the mastersheet with geo coordinates
*************************************************************
*************************************************************

use Nuts_2_Mollweide_NASA_synthesis_FINAL
spset 

merge 1:1 id using key_variables
drop _merge
save "ready_data"

drop if _ID==.
spmatrix create contiguity W, replace 
spmatrix create idistance M, replace

**************Main regressions
***Forest loss regressions
est clear 
regress pr_logging_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
spregress pr_logging_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, ml errorlag(M) vce(robust)
eststo

regress pr_logging_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
spregress pr_logging_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, ml errorlag(M) vce(robust)
eststo

regress pr_fireloss_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo 
spregress pr_fireloss_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, ml errorlag(M) vce(robust)
eststo

regress pr_fireloss_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
spregress pr_fireloss_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, ml errorlag(M) vce(robust)
eststo

*Table 2 
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest loss in 1989-2000 vs. 2001-2012) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 


****Forest gain regressions
est clear
regress pr_croppastgain1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust 
eststo
spregress pr_croppastgain1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, ml errorlag(M) vce(robust)
eststo

regress pr_croppastgain2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust 
eststo
spregress pr_croppastgain2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, ml errorlag(M) vce(robust) 
eststo 
*Table 3
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest gain in 1989-2000 vs. 2001-2012) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 


********************************************Robustness checks*************************************************************
**********Regressions dropping 1989 and 1990 as well as 2011 and 2012 
*Logging loss 
gen smaller_loggingloss_1990s =  otherloss_1991 + otherloss_1992	+ otherloss_1993 + otherloss_1994 + otherloss_1995 + otherloss_1996	+ otherloss_1997 + otherloss_1998 + otherloss_1999 + otherloss_2000
gen smaller_loggingloss_2000s = otherloss_2001 + otherloss_2002 + otherloss_2003 + otherloss_2004	+ otherloss_2005 + otherloss_2006 + otherloss_2007 + otherloss_2008	+ otherloss_2009 + otherloss_2010

gen smaller_pr_logging_1990s = (smaller_loggingloss_1990s/totalforest85)*100
gen smaller_pr_logging_2000s = (smaller_loggingloss_2000s/totalforest00)*100

*Fire loss 
gen smaller_fireloss_1990s = fireloss_1991 + fireloss_1992 + fireloss_1993 +	fireloss_1994 + fireloss_1995 + fireloss_1996 + fireloss_1997 + fireloss_1998 + fireloss_1999 + fireloss_2000
gen smaller_fireloss_2000s = fireloss_2001	+ fireloss_2002 + fireloss_2003	+ fireloss_2004	+ fireloss_2005 +	fireloss_2006	+ fireloss_2007	+ fireloss_2008	+ fireloss_2009	+ fireloss_2010

gen smaller_pr_fireloss_1990s = (smaller_fireloss_1990s/totalforest85)*100
gen smaller_pr_fireloss_2000s = (smaller_fireloss_2000s/totalforest00)*100

*Gain over cropland or pasture
gen smaller_pr_croppastgain1990s =(((aff_2000_crop_pas/16)*10)/total_nonforest85)*100
gen smaller_pr_croppastgain2000s = ((((aff_2012_crop_past)/12)*10)/total_nonforest00)*100

est clear 
regress smaller_pr_logging_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress smaller_pr_logging_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
regress smaller_pr_fireloss_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress smaller_pr_fireloss_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest loss in 1989-2000 and 2001-2012î dropping 1989/90 and 2011/12) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

est clear 
regress smaller_pr_croppastgain1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress smaller_pr_croppastgain2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest gain in 1989-2000 and 2001-2012î dropping 1989/90 and 2011/12) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

**********Regressions excluding fires at the end of the 2000s**********
*Drop 2010 fires 
gen fireloss_2000s_wo2010 = fireloss_2001 + fireloss_2002 + fireloss_2003	+ fireloss_2004	+ fireloss_2005 +	fireloss_2006	+ fireloss_2007	+ fireloss_2008	+ fireloss_2009 + fireloss_2011 + fireloss_2012
gen pr_fireloss_2000s_wo2010 = (fireloss_2000s_wo2010/totalforest00)*100

*Drop 2010 and 2011 fires
gen fireloss_2000s_wo2010_2011 = fireloss_2001 + fireloss_2002 + fireloss_2003	+ fireloss_2004	+ fireloss_2005 +	fireloss_2006	+ fireloss_2007	+ fireloss_2008	+ fireloss_2009 + fireloss_2012
gen pr_fireloss_2000s_wo2010_2011 = (fireloss_2000s_wo2010_2011/totalforest00)*100

*Drop 2009, 2010 and 2011 fires
gen fireloss_2000s_wo09_10_11 = fireloss_2001 + fireloss_2002 + fireloss_2003	+ fireloss_2004	+ fireloss_2005 +	fireloss_2006	+ fireloss_2007	+ fireloss_2008 + fireloss_2012
gen pr_fireloss_2000s_wo09_10_11 = (fireloss_2000s_wo09_10_11/totalforest00)*100

est clear 
regress pr_fireloss_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
regress pr_fireloss_2000s_wo2010 fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
regress pr_fireloss_2000s_wo2010_2011 fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
regress pr_fireloss_2000s_wo09_10_11 fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest fire loss in 1989-2000 and 2001-2012î (w/o 2010; 2010/11; 2009-11)) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

**********Regressions using alternative measurement for the 2000s forest cover**********
est clear 
regress alt_pr_logging_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 alt_ln_off2001 alt_ln_totalfor00 governance autonom, robust
eststo
regress alt_pr_fireloss_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 alt_ln_off2001 alt_ln_totalfor00 governance autonom, robust
eststo
regress alt_pr_croppastgain2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 alt_ln_off2001 alt_ln_totalfor00 governance autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest loss and gain in 1989-2000 and 2001-2012î (2000s normalized over alt forest figure)) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

***********Regressions using lower threshold for fertile area (soil index <=26)
gen larger_fertile_area = ((si26_40 + si41_55 + si56_70 + si71_85 + si86_100)/(all_classes - water))*100

est clear 
regress pr_logging_1990s larger_fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress pr_logging_2000s larger_fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
regress pr_fireloss_1990s larger_fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress pr_fireloss_2000s larger_fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest loss in 1989-2000 and 2001-2012î with FA: SI $\geq$ 26) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

est clear 
regress pr_croppastgain1990s larger_fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress pr_croppastgain2000s larger_fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest gain in 1989-2000 and 2001-2012î FA: SI $\geq$ 26) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

***********Regressions using higher threshold for fertile area (soil index <=55)
gen smaller_fertile_area = ((si56_70 + si71_85 + si86_100)/(all_classes - water))*100
est clear 
regress pr_logging_1990s smaller_fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress pr_logging_2000s smaller_fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
regress pr_fireloss_1990s smaller_fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress pr_fireloss_2000s smaller_fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest loss in 1989-2000 and 2001-2012î FA: SI $\geq$ 55) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

est clear 
regress pr_croppastgain1990s smaller_fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress pr_croppastgain2000s smaller_fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest gain in 1989-2000 and 2001-2012î FA: SI $\geq$ 55) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

*******Regressions using 1992 forest officials data; original figures for total forest in the 2000 
est clear 
regress pr_logging_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1992 ln_totalfor85 governance_91_01 autonom, robust
eststo
regress pr_fireloss_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1992 ln_totalfor85 governance_91_01 autonom, robust
eststo 
regress pr_croppastgain1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1992 ln_totalfor85 governance_91_01 autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest loss and gain in 1989-2000 and 2001-2012î using 1992 for/off data) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

*******Regressions with governance scores squared 
est clear 
regress pr_logging_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 gov_sq_90 autonom, robust
eststo
regress pr_logging_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance gov_sq_00 autonom, robust
eststo
regress pr_fireloss_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 gov_sq_90 autonom, robust
eststo
regress pr_fireloss_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance gov_sq_00 autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest loss in 1989-2000 and 2001-2012î gov.score$^2$) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

est clear 
regress pr_croppastgain1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 gov_sq_90 autonom, robust
eststo
regress pr_croppastgain2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance gov_sq_00 autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest gain in 1989-2000 and 2001-2012î gov.score$^2$) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

*******Regressions for the 2000s with regional income pc adjusted to CPI 
est clear 
regress pr_logging_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_adj_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo 
regress pr_fireloss_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_adj_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
regress pr_croppastgain2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_adj_lst2001 ln_off2001 ln_totalfor00 governance autonom, robust
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Forest loss and gain in the 2000s with adjusted income data) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 


*******Regressions for models with spatial autoregressive errors using continguity matrix 
est clear 
spregress pr_logging_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, ml errorlag(W) vce(robust)
eststo
spregress pr_logging_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, ml errorlag(W) vce(robust)
eststo
spregress pr_fireloss_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, ml errorlag(W) vce(robust)
eststo
spregress pr_fireloss_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, ml errorlag(W) vce(robust)
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest loss in 1989-2000 vs. 2001-2012) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 


est clear
spregress pr_croppastgain1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom, ml errorlag(W) vce(robust)
eststo
spregress pr_croppastgain2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom, ml errorlag(W) vce(robust) 
eststo 
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest gain in 1989-2000 vs. 2001-2012) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

******Hierarchical linear models 
est clear 
mixed pr_logging_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom ///
 || zone: 
eststo
mixed pr_logging_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom ///
 || zone: 
eststo
mixed pr_fireloss_1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom ///
 || zone: 
eststo 
mixed pr_fireloss_2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom ///
 || zone: 
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest loss in 1989-2000 and 2001-2012î) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

****Forest gain regressions
est clear 
mixed pr_croppastgain1990s fertile_area access pr_evergreen_for rescaled_dist pr_upop1991 ln_lst1990 ln_off1988 ln_totalfor85 governance_91_01 autonom ///
 || zone: 
eststo
mixed pr_croppastgain2000s fertile_area access pr_evergreen_for rescaled_dist pr_upop2001 ln_lst2001 ln_off2001 ln_totalfor00 governance autonom ///
 || zone: 
eststo
esttab using auto.tex, se b(%9.3f) r2(%9.3f) ar2(%9.3f) starlevels(* .10 ** .05 *** .001) nodepvars nomtitles title(Correlates of forest gain in 1989-2000 and 2001-2012î) addnotes(ìNote: Heteroskedascity-robust standard errorsî) replace 

****Figure 1
graph twoway (lfit ln_off2001 ln_totalfor00) (scatter ln_off2001 ln_totalfor00), ylabel(,nogrid) ytitle("Log forest bureaucrats per hectare, 2001") xtitle("Log total forest area, 2000") legend(off)
graph twoway (lfit ln_off1988 ln_totalfor85) (scatter ln_off1988 ln_totalfor85), ylabel(,nogrid) ytitle("Log forest bureaucrats per hectare, 1988") xtitle("Log total forest area, 1985") legend(off)

